load IntnlData MV;
MV(MV<=0) = NaN;
NumPort = 10; MinStocks = 50;

load MSCI Return36_Monthly RetVol12_Monthly RetVol36_Monthly;
RMKT36 = Return36_Monthly(230:end-4,2:4);
Vol12 = RetVol12_Monthly(230:end-4,2:4);
Vol36 = RetVol36_Monthly(230:end-4,2:4);
UP = false(T,55); DOWN = false(T,55); HIVOL = false(T,55); LOVOL = false(T,55);
UP(:,[55 50 51]) = RMKT36>=0;  DOWN(:,[55 50 51]) = RMKT36<0; 
HIVOL(:,[55 50 51]) = Vol12 > Vol36;  LOVOL(:,[55 50 51]) = Vol12 <= Vol36; 
clear Return36_Monthly RetVol12_Monthly RetVol36_Monthly RMKT36 Vol12 Vol36 tmp;

TFirst = find(ym==199303); THalf = find(ym==200612);
for cntry = [55 50 51]
    RPort = NaN(NumPort+1,2,T);
    NPort = NaN(NumPort+1,T);
    for t = TFirst-2:T-2
        R11_t = Return11(:,t);
        MV_t1 = MV(:,t+1);
        idx0 = isfinite(R11_t) & isfinite(MV_t1) & CountryRegionDummy(:,cntry) & ~IsREIT & ~IsMicro(:,t,cntry);
        if sum(idx0) < MinStocks*NumPort %number of stocks at this stage should be enough to populate all portfolios
            continue;
        end

        R11_t(idx0) = winsorize(R11_t(idx0),WinsorLevelX);
        RNext = Return(:,t+2);
        if WinsorizeYRet
            RNext(idx0) = winsorize(RNext(idx0),WinsorLevelY);
        end
        MVNext = MV(:,t+1);    MVNext(idx0) = winsorize(MVNext(idx0),WinsorLevelY);

        % neutralize country effects
        for cc = 1:49
            idx_sort_cc = CountryRegionDummy(:,cc) & idx0;
            if sum(idx_sort_cc) > 0
                R11_t(idx_sort_cc) = R11_t(idx_sort_cc) - mean(R11_t(idx_sort_cc),'omitnan');
            end
        end

        Rprev_sort = sort(R11_t(idx0));
        NinPort = length(Rprev_sort)/NumPort;
        BrkPts_n = Rprev_sort([1 floor(NinPort*(1:NumPort))]); BrkPts_n(end) = Inf;

        for n = 1:NumPort
            idx_n = R11_t>=BrkPts_n(n) & R11_t<BrkPts_n(n+1) & idx0;
            if sum(idx_n) >= MinStocks %ensure that number of stocks in this portfolio is enough
                ret = RNext(idx_n); mv = MVNext(idx_n);
                idx = isfinite(ret) & isfinite(mv);
                ret = ret(idx); mv = mv(idx);
                if sum(idx) >= MinStocks % ensure again that we have returns for enough stocks
                    RPort(n,1,t+2) = sum(ret.*mv)/sum(mv);
                    RPort(n,2,t+2) = mean(ret);
                    NPort(n,t+2) = sum(idx);
                end
            end
        end
    end
    % take data only where corner portfolios are populated
    tmp = [RPort(1,1,:) RPort(NumPort,1,:)];
    idx = all(isfinite(tmp),2);
    RPort(:,:,~idx) = NaN; NPort(:,~idx) = NaN;
    RPort(NumPort+1,:,:) = RPort(NumPort,:,:) - RPort(1,:,:);

    ToDisp1 = NaN(6,14); ToDisp2 = NaN(6,14);
    ewvw = 1;

    for sample = 1:3
        switch sample
            case 1, T1 = TFirst; T2 = T;
            case 2, T1 = TFirst; T2 = THalf;
            case 3, T1 = THalf+1; T2 = T;
        end
        for pp = 1:3 % Loser, Winner, Winner-Loser
            if pp==1; Y = squeeze(RPort(1,ewvw,T1:T2)); elseif pp==2; Y = squeeze(RPort(NumPort,ewvw,T1:T2)); else Y = squeeze(RPort(end,ewvw,T1:T2)); end

            X = [UP(T1:T2,cntry) DOWN(T1:T2,cntry)];
            x = X(1:end-1,:);
            y = Y(2:end);

            idx = all(isfinite(x),2) & isfinite(y);
            y = y(idx); x = x(idx,:);
            [B,Tols] = myols_inside(x,y);
            colnum = (sample-1)*5 + 1;
            ToDisp1(1,colnum) = sum(x(:,1));
            ToDisp1(3,colnum) = sum(x(:,2));

            colnum = (sample-1)*5 + 1 + pp;
            ToDisp1(1,colnum) = 1200*B(1); ToDisp1(2,colnum) = Tols(1);
            ToDisp1(3,colnum) = 1200*B(2); ToDisp1(4,colnum) = Tols(2);
            ToDisp1(5,colnum) = 1200*B(3); ToDisp1(6,colnum) = Tols(3);

            X = [HIVOL(T1:T2,cntry) LOVOL(T1:T2,cntry)];
            x = X(1:end-1,:);
            y = Y(2:end);

            idx = all(isfinite(x),2) & isfinite(y);
            y = y(idx); x = x(idx,:);
            [B,Tols] = myols_inside(x,y);
            colnum = (sample-1)*5 + 1;
            ToDisp2(1,colnum) = sum(x(:,1));
            ToDisp1(3,colnum) = sum(x(:,2));

            colnum = (sample-1)*5 + 1 + pp;
            ToDisp2(1,colnum) = 1200*B(1); ToDisp2(2,colnum) = Tols(1);
            ToDisp2(3,colnum) = 1200*B(2); ToDisp2(4,colnum) = Tols(2);
            ToDisp2(5,colnum) = 1200*B(3); ToDisp2(6,colnum) = Tols(3);
        end
    end

    D = [ToDisp1 NaN(6,1) ToDisp2];
    if cntry == 55
        writematrix(D,FileName, 'Sheet', SheetName, 'Range','b5:ad10', 'AutoFitWidth',false);
    elseif cntry == 50
        writematrix(D,FileName, 'Sheet', SheetName, 'Range','b12:ad17', 'AutoFitWidth',false);
    elseif cntry == 51
        writematrix(D,FileName, 'Sheet', SheetName, 'Range','b19:ad24', 'AutoFitWidth',false);
    end
end
  
%---------------------------------------------------------------------------
function [B,Tols] = myols_inside(X,Y)
idx = all(isfinite(X),2) & isfinite(Y);
X = X(idx,:); Y = Y(idx);

[T,K] = size(X);
B = X\Y; E = Y-X*B;
Q = inv(X'*X);
s = E'*E/(T-K);
V_ols = s*Q;

w = [1 0; 0 1; 1 -1];
B = w*B;
V_ols = w*V_ols*w';

Tols = B./sqrt(diag(V_ols));
end
